Dicom loading and visualization¶

Cargar Archivos DICOM¶

Cargo los archivos DICOM desde el directorio especificado. Trabajo inicialmente con un subconjunto de los archivos para testear.

¿Con que CT trabajo? No me ha quedado claro qué subconjunto de archivos usar!

  1. ct_directory = "C:/Users/Pedro/Documents/GitHub/Medial_Image_Pydicom/HCC_007/manifest-1643035385102/HCC-TACE-Seg/HCC_007/03-14-1998-NA-CT ABDOMEN WWO CONT-48924/4.000000-Recon 2 3 PHASE LIVER ABD-87008"
  2. ct_directory_2 = "C:/Users/Pedro/Documents/GitHub/Medial_Image_Pydicom/HCC_007/manifest-1643035385102/HCC-TACE-Seg/HCC_007/03-14-1998-NA-CT ABDOMEN WWO CONT-48924/2.000000-PRE LIVER-48636"

Libraries¶

In [ ]:
import os
import pydicom
import numpy as np
import matplotlib.pyplot as plt
from typing import List
import cv2
import matplotlib.colors
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation

Aqui alojamos las funciones¶

In [ ]:
def load_dicom_files(directory: str, file_extension: str = '.dcm'):
    """Load DICOM files from a specified directory that match the given file extension."""
    return [pydicom.dcmread(os.path.join(directory, f)) for f in os.listdir(directory) if f.endswith(file_extension)]
    
def sort_dicom_files_by_instance_number(dicom_files: List[pydicom.Dataset]) -> List[pydicom.Dataset]:
    """Sort DICOM files by the Instance Number."""
    return sorted(dicom_files, key=lambda x: int(x.InstanceNumber))

def sort_dicom_files_by_position(dicom_files):
    """Ordena archivos DICOM según la posición del paciente en la imagen."""
    return sorted(dicom_files, key=lambda x: float(x.ImagePositionPatient[2]))

def sort_dicom_files(dicom_files):
    """Sort DICOM files by available headers. Handles both CT and segmentation files."""
    def get_sort_key(file):
        # Define sorting keys: AcquisitionNumber if available, otherwise ImagePositionPatient
        acquisition_number = getattr(file, 'AcquisitionNumber', None)
        position_patient = getattr(file, 'ImagePositionPatient', None)

        # Convert ImagePositionPatient to a sortable form if available
        if position_patient:
            position_patient = float(position_patient[2])
        return (int(acquisition_number) if acquisition_number else 0, position_patient if position_patient else 0)
    '''
    Primero ordena por AcquisitionNumber y, dentro de cada número de adquisición, 
    ordena por la posición del paciente a lo largo del eje Z (ImagePositionPatient[2]).'''
    return sorted(dicom_files, key=get_sort_key)

def check_acquisition_uniformity(dicom_files: List[pydicom.Dataset]):
    """Check if all DICOM files have the same 'Acquisition Number'."""
    acquisition_numbers = {img.AcquisitionNumber for img in dicom_files}
    if len(acquisition_numbers) > 1:
        print("Multiple acquisition numbers found. Ensure all CT images are from a single acquisition.")
    
def get_pixel_arrays(dicom_files: List[pydicom.Dataset]) -> np.ndarray:
    """Extract pixel arrays from sorted DICOM files and stack them into a single numpy array."""
    return np.stack([file.pixel_array for file in dicom_files])

def display_images(ct_image, seg_image):
    """Display CT and Segmentation images side-by-side."""
    plt.figure(figsize=(12, 6))

    plt.subplot(1, 2, 1)
    plt.imshow(ct_image, cmap='gray', aspect='auto')
    plt.title('CT Slice')

    plt.subplot(1, 2, 2)
    plt.imshow(seg_image, cmap='gray', aspect='auto')
    plt.title('Segmentation Slice')
    plt.show()

def create_mip(ct_volume: np.ndarray):
    """Create and visualize Maximum Intensity Projection (MIP) of a CT volume."""
    mip_axial = np.max(ct_volume, axis=0)
    mip_coronal = np.max(ct_volume, axis=1)
    mip_sagittal = np.max(ct_volume, axis=2)
    
    projections = [mip_axial, mip_coronal, mip_sagittal]
    titles = ['MIP Axial', 'MIP Coronal', 'MIP Sagittal']

    plt.figure(figsize=(15, 5))
    for i, mip in enumerate(projections):
        plt.subplot(1, 3, i + 1)
        vmin, vmax = np.percentile(mip, [2, 98])
        plt.imshow(mip, cmap='gray', aspect='auto', vmin=vmin, vmax=vmax)
        plt.title(titles[i])
        plt.axis('off')
    plt.show()

def visualize_planes(ct_volume: np.ndarray):
    """Visualize axial, coronal, and sagittal planes of a CT volume."""
    coronal_plane = ct_volume[:, :, ct_volume.shape[2] // 2]
    sagittal_plane = ct_volume[:, ct_volume.shape[1] // 2, :]
    axial_plane = ct_volume[ct_volume.shape[0] // 2, :, :]
    
    plt.figure(figsize=(15, 5))
    for i, (plane, title) in enumerate(zip([axial_plane, coronal_plane, sagittal_plane],
                                           ['Axial Plane', 'Coronal Plane', 'Sagittal Plane'])):
        plt.subplot(1, 3, i+1)
        plt.imshow(plane, cmap='gray', aspect='auto')
        plt.title(title)
        plt.axis('off')
    plt.show()

def apply_window(image: np.ndarray, window_center: int, window_width: int) -> np.ndarray:
    """Apply window leveling to an image based on the given center and width."""
    lower_bound = window_center - 0.5 * window_width
    upper_bound = window_center + 0.5 * window_width
    return np.clip((image - lower_bound) / window_width, 0, 1)
In [ ]:
# Data paths
ct_directory = "C:/Users/Pedro/Documents/GitHub/Medial_Image_Pydicom/HCC_007/manifest-1643035385102/HCC-TACE-Seg/HCC_007/03-14-1998-NA-CT ABDOMEN WWO CONT-48924/4.000000-Recon 2 3 PHASE LIVER ABD-87008"
seg_directory = "C:/Users/Pedro/Documents/GitHub/Medial_Image_Pydicom/HCC_007/manifest-1643035385102/HCC-TACE-Seg/HCC_007/12-27-1997-NA-AP LIVER PRT WWO-67834/300.000000-Segmentation-39839/1-1.dcm"
#ct_directory_2 = "C:/Users/Pedro/Documents/GitHub/Medial_Image_Pydicom/HCC_007/manifest-1643035385102/HCC-TACE-Seg/HCC_007/03-14-1998-NA-CT ABDOMEN WWO CONT-48924/2.000000-PRE LIVER-48636"

Cargamos las imagenes y las ordenamos¶

Cargamos y ordenamos las imagenes CT

In [ ]:
# Load and process DICOM CT images
ct_files = load_dicom_files(ct_directory)
sorted_ct_files = sort_dicom_files(ct_files)

Cargamos la segmentation image

In [ ]:
# Load and process segmentation image
segmentation_image = pydicom.dcmread(seg_directory).pixel_array
sorted_seg_img = sort_dicom_files(segmentation_image)
mid_slice_index = segmentation_image.shape[0] // 2

Simple visualización para saber que todo ha ido bien

In [ ]:
# Visualize a slice of CT and Segmentation images
ct_images = get_pixel_arrays(sorted_ct_files)
display_images(ct_images[8], segmentation_image[8])
No description has been provided for this image

Estudio de encabezados¶

In [ ]:
# Print details and dimensions
print("Details of CT images sorted by Image Position Patient:")
for img in sorted_ct_files:
    print(f"Acquisition Number: {img.AcquisitionNumber}, Image Position Patient: {img.ImagePositionPatient}")

print("Dimensiones del arreglo de segmentación:", segmentation_image.shape)
print("Dimensiones del arreglo de imágenes CT:", ct_images.shape)
Details of CT images sorted by Image Position Patient:
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -260.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -257.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -255.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -252.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -250.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -247.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -245.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -242.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -240.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -237.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -235.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -232.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -230.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -227.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -225.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -222.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -220.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -217.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -215.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -212.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -210.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -207.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -205.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -202.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -200.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -197.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -195.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -192.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -190.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -187.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -185.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -182.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -180.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -177.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -175.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -172.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -170.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -167.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -165.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -162.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -160.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -157.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -155.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -152.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -150.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -147.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -145.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -142.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -140.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -137.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -135.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -132.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -130.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -127.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -125.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -122.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -120.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -117.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -115.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -112.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -110.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -107.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -105.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -102.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -100.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -97.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -95.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -92.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -90.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -87.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -85.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -82.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -80.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -77.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -75.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -72.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -70.000000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -67.500000]
Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -65.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -350.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -347.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -345.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -342.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -340.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -337.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -335.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -332.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -330.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -327.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -325.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -322.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -320.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -317.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -315.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -312.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -310.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -307.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -305.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -302.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -300.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -297.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -295.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -292.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -290.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -287.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -285.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -282.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -280.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -277.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -275.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -272.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -270.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -267.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -265.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -262.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -260.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -257.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -255.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -252.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -250.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -247.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -245.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -242.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -240.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -237.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -235.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -232.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -230.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -227.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -225.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -222.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -220.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -217.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -215.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -212.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -210.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -207.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -205.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -202.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -200.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -197.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -195.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -192.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -190.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -187.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -185.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -182.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -180.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -177.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -175.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -172.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -170.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -167.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -165.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -162.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -160.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -157.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -155.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -152.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -150.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -147.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -145.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -142.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -140.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -137.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -135.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -132.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -130.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -127.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -125.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -122.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -120.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -117.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -115.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -112.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -110.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -107.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -105.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -102.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -100.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -97.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -95.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -92.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -90.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -87.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -85.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -82.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -80.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -77.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -75.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -72.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -70.000000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -67.500000]
Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -65.000000]
Dimensiones del arreglo de segmentación: (300, 512, 512)
Dimensiones del arreglo de imágenes CT: (194, 512, 512)

Confirmamos que hay más de una adquisición

Que hago aqui? Hay más de una adquisición pero solo debe haber una. En los obejtivos pone "Asegurar que solo haya una adquisicion".

In [ ]:
check_acquisition_uniformity(sorted_ct_files)
Multiple acquisition numbers found. Ensure all CT images are from a single acquisition.

Aplicamos windowing a las imagenes para mejorarlas¶

In [ ]:
# Windowing parameters
window_center, window_width = 40, 400
windowed_ct_images = np.array([apply_window(img, window_center, window_width) for img  in ct_images])
In [ ]:
# Visualize different planes
visualize_planes(windowed_ct_images)
No description has been provided for this image
In [ ]:
# Visualize images and segmentations
for i in range(10):  # Example for the first 10 slices
    display_images(windowed_ct_images[i], segmentation_image[i])
print("Segmentation array loaded with shape:", segmentation_image.shape)
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
Segmentation array loaded with shape: (300, 512, 512)

ROI¶

In [ ]:
def apply_threshold(ct_slice, threshold=0.1):
    """Aplicar umbralización para separar el ROI del fondo."""
    _, thresh = cv2.threshold(ct_slice, threshold, 255, cv2.THRESH_BINARY)
    return thresh
'''
PREGUNTA: Tenemos que excluir aquellas zonasa que son demasiado pequeñas?
 Se podría implementar un filtro por tamaño,
 conservando solo las ROIs que tengan un número mínimo de píxeles, 
 lo que indica que son regiones significativas.
 
 Ej. unique, counts = np.unique(labels_im, return_counts=True)
     filtered_labels = np.zeros_like(labels_im)
     for label, count in zip(unique, counts):
        if count > min_pixel_count:  # Filtrar ROIs pequeñas
            filtered_labels[labels_im == label] = label
'''
def find_rois(thresh_image): 
    """Encontrar y etiquetar ROIs en la imagen umbralizada."""
    num_labels, labels_im = cv2.connectedComponents(thresh_image.astype(np.uint8))
    return num_labels, labels_im

def visualize_rois(ct_slice, labels_im):
    """Visualizar las ROIs sobre la imagen original."""
    plt.figure(figsize=(10, 5))
    plt.subplot(1, 2, 1)
    plt.imshow(ct_slice, cmap='gray')
    plt.title('Original CT Slice')
    plt.subplot(1, 2, 2)
    plt.imshow(labels_im, cmap='nipy_spectral')
    plt.title('Detected ROIs')
    plt.show()

ct_slice = windowed_ct_images[0]

# Análisis de datos de la imagen
print("Estadísticas de la imagen:")
print(f"Valor mínimo: {np.min(ct_slice)}, Valor máximo: {np.max(ct_slice)}")
print(f"Media: {np.mean(ct_slice)}, Desviación estándar: {np.std(ct_slice)}")

# Procesamiento
threshold = 0.43
thresh = apply_threshold(ct_slice, threshold)
print(f"Porcentaje de píxeles sobre el umbral: {np.sum(thresh > 0) / thresh.size * 100:.2f}%")

num_labels, labels_im = find_rois(thresh)
print(f"Número de ROIs detectadas: {num_labels - 1}")  # Excluyendo el fondo

# Visualización
visualize_rois(ct_slice, labels_im)

# Análisis de componentes conectados
'''
unique, counts = np.unique(labels_im, return_counts=True)
for roi_label, count in zip(unique, counts):
    if roi_label != 0:  # Excluyendo el fondo
        print(f"ROI Label: {roi_label}, Pixel Count: {count}")
'''
Estadísticas de la imagen:
Valor mínimo: 0.0, Valor máximo: 1.0
Media: 0.2165085220336914, Desviación estándar: 0.2696629915027511
Porcentaje de píxeles sobre el umbral: 25.47%
Número de ROIs detectadas: 208
No description has been provided for this image
Out[ ]:
'\nunique, counts = np.unique(labels_im, return_counts=True)\nfor roi_label, count in zip(unique, counts):\n    if roi_label != 0:  # Excluyendo el fondo\n        print(f"ROI Label: {roi_label}, Pixel Count: {count}")\n'

MIP¶

In [ ]:
def verify_alignment(ct_dcm, seg_dcm):
    """Verify if the DICOM files have matching slice locations."""
    if hasattr(ct_dcm, 'SliceLocation') and hasattr(seg_dcm, 'SliceLocation'):
        if abs(ct_dcm.SliceLocation - seg_dcm.SliceLocation) > 0.1:  # Assuming a tolerance of 0.1
            print(f"Mismatch: CT {ct_dcm.SliceLocation} vs. Seg {seg_dcm.SliceLocation}")
            return False
        print("CT and Segmentation are aligned.")
        return True
    print("One or both files do not contain SliceLocation data for comparison.")
    return False

def visualize_ct_with_segmentation(ct_slice, seg_slice):
    """Visualize CT slice with segmentation overlay."""
    fig, ax = plt.subplots(figsize=(10, 8))
    ax.imshow(ct_slice, cmap='gray')
    colored_cmap = create_transparent_cmap(plt.cm.viridis, 0.3)
    ax.imshow(seg_slice, cmap=colored_cmap, interpolation='none')
    ax.set_title('CT Slice with Segmentation Overlay')
    plt.axis('off')
    plt.show()

def create_transparent_cmap(cmap, transparency):
    """Create a colormap with transparency for overlay."""
    colors = cmap(np.linspace(0, 1, 256))
    colors[:, -1] = np.linspace(0, transparency, 256)
    return matplotlib.colors.ListedColormap(colors)

def calculate_mip(ct_volume, axis):
    """Calculate the Maximum Intensity Projection (MIP) along a given axis."""
    return np.max(ct_volume, axis=axis)

def calculate_and_display_mips(ct_volume):
    """Calculate MIPs for all three planes and display them with a lilac color theme."""
    mip_axial = calculate_mip(ct_volume, axis=0)
    mip_coronal = calculate_mip(ct_volume, axis=1)
    mip_sagittal = calculate_mip(ct_volume, axis=2)
    display_mip_images(mip_axial, mip_coronal, mip_sagittal)

def display_mip_images(mip_axial, mip_coronal, mip_sagittal):
    """Display MIP images for axial, coronal, and sagittal planes with customized color."""
    plt.figure(figsize=(18, 6))

    # Colormap personalizado que tiende al lila
    cmap = plt.get_cmap('twilight')

    # Definir los percentiles para ajustar los valores de visualización de cada MIP
    vmin_a, vmax_a = np.percentile(mip_axial, [2, 98])
    vmin_c, vmax_c = np.percentile(mip_coronal, [2, 98])
    vmin_s, vmax_s = np.percentile(mip_sagittal, [2, 98])

    # Axial
    plt.subplot(1, 3, 1)
    plt.imshow(mip_axial, cmap=cmap, aspect='auto', vmin=vmin_a, vmax=vmax_a)
    plt.title('MIP Axial')
    plt.axis('off')

    # Coronal
    plt.subplot(1, 3, 2)
    plt.imshow(mip_coronal, cmap=cmap, aspect='auto', vmin=vmin_c, vmax=vmax_c)
    plt.title('MIP Coronal')
    plt.axis('off')

    # Sagittal
    plt.subplot(1, 3, 3)
    plt.imshow(mip_sagittal, cmap=cmap, aspect='auto', vmin=vmin_s, vmax=vmax_s)
    plt.title('MIP Sagittal')
    plt.axis('off')

    plt.tight_layout()
    plt.show()
In [ ]:
# Example usage within your data processing pipeline:
'''for ct_dcm, seg_dcm in zip(ct_images, segmentation_image):
    if verify_alignment(ct_dcm, seg_dcm):
        visualize_ct_with_segmentation(ct_dcm.pixel_array, seg_dcm.pixel_array)'''

calculate_and_display_mips(ct_images)
No description has been provided for this image

3D¶

In [ ]:
def plot_3d(image, threshold=-300, rotation_angle=360, interval=50):
    """Plot 3D visualization of CT data."""
    p = image.transpose(2, 1, 0)
    p = p.clip(min=threshold)
    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(111, projection='3d')
    pos = np.where(p > threshold)
    scat = ax.scatter(pos[0], pos[1], pos[2], c=p[pos], cmap='viridis', s=1)
    ax.view_init(30, 185)
    ani = FuncAnimation(fig, lambda frame: ax.view_init(30, frame), frames=np.arange(0, rotation_angle, 2), interval=interval)
    plt.colorbar(scat, ax=ax, orientation='vertical')
    plt.show()
    return ani

plot_3d(ct_images, threshold=100)  # Adjust threshold as needed
No description has been provided for this image
Out[ ]:
<matplotlib.animation.FuncAnimation at 0x15f70622760>

Pregunta: Como quito la linea de atras? Esa columna que se ve cerca del eje y

Prueba: Superposición de las CT sobre la segmentada¶

In [ ]:
def visualize_ct_with_segmentation(ct_slice, seg_slice):
    
    # Crear figura y ejes
    fig, ax = plt.subplots(figsize=(10, 8))
    ax.imshow(ct_slice, cmap='gray')  # Visualizar la imagen de CT en escala de grises
    
    # Configuración del colormap para la segmentación con transparencia
    cmap = plt.cm.viridis
    colors = cmap(np.linspace(0, 1, 256))
    colors[:, -1] = np.linspace(0, 0.3, 256)  # Configuración de la transparencia
    colored_cmap = matplotlib.colors.ListedColormap(colors)

    # Superponer la segmentación en la imagen de CT
    seg_im = ax.imshow(seg_slice, cmap=colored_cmap, interpolation='none')
    
    ax.set_title('CT Slice with Segmentation Overlay')
    plt.axis('off')
    plt.show()

# Iterar sobre cada slice en el conjunto de datos de CT y Segmentación
for i in range(5):
    ct_slice = ct_images[i+5]
    seg_slice = segmentation_image[i+5]
    visualize_ct_with_segmentation(ct_slice, seg_slice)
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

Image Coregistration¶

Cargamos las Imagenes DICOM¶

In [ ]:
reference_image_path = "C:/Users/Pedro/Documents/GitHub/Medial_Image_Pydicom/Project/icbm_avg_152_t1_tal_nlin_symmetric_VI.dcm"
input_images_path = "C:/Users/Pedro/Documents/GitHub/Medial_Image_Pydicom/Project/RM_Brain_3D-SPGR"
In [ ]:
input_images = sort_dicom_files_by_position(load_dicom_files(input_images_path))
reference_image = pydicom.dcmread(reference_image_path)
In [ ]:
def multiply_quaternions(
        q1: tuple[float, float, float, float],
        q2: tuple[float, float, float, float]
        ) -> tuple[float, float, float, float]:
    """ Multiply two quaternions, expressed as (1, i, j, k). """
    return (
        q1[0] * q2[0] - q1[1] * q2[1] - q1[2] * q2[2] - q1[3] * q2[3],
        q1[0] * q2[1] + q1[1] * q2[0] + q1[2] * q2[3] - q1[3] * q2[2],
        q1[0] * q2[2] - q1[1] * q2[3] + q1[2] * q2[0] + q1[3] * q2[1],
        q1[0] * q2[3] + q1[1] * q2[2] - q1[2] * q2[1] + q1[3] * q2[0]
    )


def conjugate_quaternion(
        q: tuple[float, float, float, float]
        ) -> tuple[float, float, float, float]:
    """ Multiply two quaternions, expressed as (1, i, j, k). """
    return (
        q[0], -q[1], -q[2], -q[3]
    )

def mutual_information(hgram):
    """ Calculate the mutual information between two images """
    pxy = hgram / float(np.sum(hgram))
    px = np.sum(pxy, axis=1)
    py = np.sum(pxy, axis=0)
    px_py = px[:, None] * py[None, :]
    nzs = pxy > 0
    return np.sum(pxy[nzs] * np.log(pxy[nzs] / px_py[nzs]))


def translation(
        point: tuple[float, float, float],
        translation_vector: tuple[float, float, float]
        ) -> tuple[float, float, float]:
    """ Perform translation of `point` by `translation_vector`. """
    x, y, z = point
    v1, v2, v3 = translation_vector
    return (x+v1, y+v2, z+v3)

def transform_point(point, t, q):
    """Apply translation and rotation to a point using quaternions."""
    # Translate
    point = np.array(point) + np.array(t)
    # Rotate using quaternion
    p = (0,) + tuple(point)
    p_rotated = multiply_quaternions(multiply_quaternions(q, p), conjugate_quaternion(q))
    return p_rotated[1:]

def axial_rotation(
        point: tuple[float, float, float],
        angle_in_rads: float,
        axis_of_rotation: tuple[float, float, float]) -> tuple[float, float, float]:
    """ Perform axial rotation of `point` around `axis_of_rotation` by `angle_in_rads`. """
    x, y, z = point
    v1, v2, v3 = axis_of_rotation

    # Normalize axis of rotation to avoid restrictions on optimizer
    v_norm = math.sqrt(sum([coord ** 2 for coord in [v1, v2, v3]]))
    v1, v2, v3 = v1 / v_norm, v2 / v_norm, v3 / v_norm

    #   Quaternion associated to point.
    p = (0, x, y, z)

    #   Quaternion associated to axial rotation.
    cos, sin = math.cos(angle_in_rads / 2), math.sin(angle_in_rads / 2)
    q = (cos, sin * v1, sin * v2, sin * v3)

    #   Quaternion associated to image point
    q_star = conjugate_quaternion(q)
    p_prime = multiply_quaternions(q, multiply_quaternions(p, q_star))
    
    #   Interpret as 3D point (i.e. drop first coordinate)
    return p_prime[1], p_prime[2], p_prime[3]

def ssd(reference, transformed):
    """Compute Sum of Squared Differences between two images."""
    return np.sum((reference - transformed) ** 2)

def transform_image(image, params):
    """Apply translation and rotation to the image based on params."""
    # Unpack translation and rotation parameters
    t1, t2, t3, angle, axis = params
    translated = translation(image, (t1, t2, t3))
    rotated = axial_rotation(translated, angle, axis)
    return rotated

def objective_function(params, input_image, reference_image):
    """Objective function to minimize (SSD between reference and transformed input)."""
    transformed_image = transform_image(input_image, params)
    return ssd(reference_image, transformed_image)

def transform_image(image, params):
    """Apply translation and rotation to the image based on params."""
    # Unpack translation and rotation parameters
    t1, t2, t3, angle, axis = params
    translated = translation(image, (t1, t2, t3))
    rotated = axial_rotation(translated, angle, axis)
    return rotated

def objective_function(params, input_image, reference_image):
    """Objective function to minimize (SSD between reference and transformed input)."""
    transformed_image = transform_image(input_image, params)
    return ssd(reference_image, transformed_image)
In [ ]:
# Initial transformation parameters: translation (tx, ty, tz) and rotation (angle, axis)
initial_params = [0, 0, 0, 0, (0, 0, 1)]  # example initial guess

# Optimization to find the best parameters
result = minimize(
    objective_function, 
    initial_params, 
    args=(input_images, reference_image),
    method='BFGS'
)

# Apply the optimized transformation to the input image
optimized_params = result.x
transformed_image = transform_image(input_images, optimized_params)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[220], line 5
      2 initial_params = [0, 0, 0, 0, (0, 0, 1)]  # example initial guess
      4 # Optimization to find the best parameters
----> 5 result = minimize(
      6     objective_function, 
      7     initial_params, 
      8     args=(input_images, reference_image),
      9     method='BFGS'
     10 )
     12 # Apply the optimized transformation to the input image
     13 optimized_params = result.x

File c:\Users\Pedro\anaconda3\envs\AIV_P1\lib\site-packages\scipy\optimize\_minimize.py:530, in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
     51 def minimize(fun, x0, args=(), method=None, jac=None, hess=None,
     52              hessp=None, bounds=None, constraints=(), tol=None,
     53              callback=None, options=None):
     54     """Minimization of scalar function of one or more variables.
     55 
     56     Parameters
   (...)
    528 
    529     """
--> 530     x0 = np.atleast_1d(np.asarray(x0))
    532     if x0.ndim != 1:
    533         raise ValueError("'x0' must only have one dimension.")

ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (5,) + inhomogeneous part.
In [ ]:
from scipy.optimize import least_squares

def registration_error(params, ref_points, input_points):
    """Calculate error between reference points and transformed input points."""
    t = params[:3]
    q = (params[3],) + tuple(params[4:7])  # Quaternion parameters
    transformed_points = [transform_point(inp, t, q) for inp in input_points]
    return np.linalg.norm(np.array(transformed_points) - np.array(ref_points), axis=1)

# Puntos de referencia random
ref_landmarks = np.array([[10, 20, 30], [40, 50, 60], [70, 80, 90]])
input_landmarks = np.array([[12, 22, 32], [42, 52, 62], [72, 82, 92]])

# Initial parameters: [tx, ty, tz, qw, qx, qy, qz]
initial_params = [0, 0, 0, 1, 0, 0, 0]  # Identity quaternion and no translation

result = least_squares(registration_error, initial_params, args=(ref_landmarks, input_landmarks))

print("Optimal parameters:", result.x)
Optimal parameters: [-0.01442782 -0.01436382 -0.01429953  0.98456473  0.00100211 -0.00199494
  0.000993  ]

Visualization of Talamus Region¶